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Computational  fluid  dynamics  (CFD)  tools  for  the  modelling  of  solid  oxide  fuel  cells  (SOFC),  solid  oxide 
electrolyser  cells  (SOEC)  or  solid  oxide  regenerative  fuel  cells  (SORFC)  nearly  always  require  a  fitting  process 
prior  to  its  application  for  cell  design  or  optimisation  purposes.  In  this  fitting,  a  set  of  experimental  data 
is  used  to  guess  the  value  of  those  parameters  of  the  model  that  cannot  be  either  modelled  or  measured 
experimentally.  This  is  crucially  the  case  of  the  charge  transfer  coefficients  (aba,afa,abc,abc)  and  the 
exchange  current  densities  (i00,  iac)  in  the  Butler-Volmer  equation  (i.e.  electrochemical  model). 

The  fitting  of  the  electrochemical  parameters  in  the  SOFC,  SOEC  and  SORFC  modelling  literature  is 
reviewed  in  this  work.  It  is  found  that  this  process  is  only  vaguely  discussed,  if  mentioned  at  all.  In  the 
authors'  opinion,  this  practice  contributes  with  uncertainty  rather  than  guidance,  since  this  fitting  process  is 
of  utmost  significance  for  making  reliable  quantitative  predictions. 

In  this  work,  we  further  introduce  a  comprehensive  model  for  the  simulation  of  solid  oxide  regenerative 
fuel  cells,  i.e.  a  model  that  simulates,  without  any  ad  hoc  adjustments  or  tuning,  both  the  SOFC  and  SOEC 
modes.  We  also  describe  in  detail  how  the  electrochemical  parameters  are  fitted,  and  discuss  the  applicability 
of  the  values  commonly  used  in  the  literature  for  these  fitted  parameters  and  their  proper  validation.  Finally, 
the  validity  of  the  proposed  model  and  fitted  parameters  is  shown  by  comparison  of  the  numerical  results 
with  experimental  data. 
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1.  Introduction 

Hydrogen  has  been  identified  as  a  promising  energy  carrier:  it 
can  be  cleanly  produced  in  a  water  electrolyser  driven  by  renew¬ 
able  energy;  it  may  be  then  stored  or  transported,  and  eventually 
converted  back  into  water  and  electricity  in  a  fuel  cell.  A  Solid 
Oxide  Regenerative  Fuel  Cell  (SORFC)  is  an  electrochemical  device 
that  allows  both  operating  modes,  i.e.  it  can  operate  reversibly 
either  as  a  fuel  cell  or  as  an  electrolyser. 

An  SORFC  operating  as  a  solid  oxide  fuel  cell  (SOFC)  generates 
electricity,  heat  and  water  from  hydrogen  and  oxygen.  An  SORFC 
operating  as  a  solid  oxide  electrolyser  cell  (SOEC)  produces  hydrogen 
and  oxygen  when  water,  heat  and  electricity  are  supplied. 

The  reversibility  of  the  solid  oxide  cells  was  first  proved  in  the 
1990s  for  both  tubular  and  planar  geometries  [3-5];  however, 
research  in  this  field  subsided  due  to  the  low  prices  of  fossil  fuels. 
In  recent  years  environmental,  economic  and  geopolitical  concerns 
over  fossil  fuels  have  rekindled  the  interest  in  this  technology.  This 
is  evinced  by  the  ongoing  research  work  for  both  planar  [6-8]  and 
microtubular  solid  oxide  regenerative  fuel  cells  [9,10], 

Numerical  modelling  is  a  primary  tool  to  understand,  and 
optimise,  the  operation  of  solid  oxide  cells,  either  in  fuel-cell  or 
electrolyser  mode.  The  modelling  of  either  SOFCs  or  SOECs  has 
been  the  subject  of  many  papers  [11-19],  Often  such  models  use 
computational  fluid  dynamics  (CFD)  to  simulate  all  the  relevant 
mass,  heat  and  charge  transfer  processes.  Wang  et  al.  [1]  present 
an  overview  of  the  several  modelling  alternatives,  including: 
physical  models  (by  which  they  mean  those  that  represent, 
mathematically,  the  underlying  physics,  be  it  in  zero,  one,  multiple 
spatial  dimensions);  equivalent  circuit  models  (based  on  electro¬ 
chemical  impedance  spectroscopy,  or  EIS,  measurements);  and 
“gray-”  and  “black-box”  models  (such  as  Artificial  Neural  Networks 


and  neuro-fuzzy  systems).  The  review  paper  by  Hajimolana  et  al. 
[2]  is  a  systematic  inventory  of  the  main  submodels  proposed 
for  the  mathematical  representation  of  all  the  relevant  physical 
processes  in  the  several  spatial  domains  (gas  channels,  electrodes, 
electrolyte,  interconnects). 

References  reporting  the  modelling  of  both  SOFC  and  SOEC 
operating  modes  (i.e.  SORFC)  with  a  single  CFD  model  that  works 
satisfactorily  in  both  situations  are  very  scarce.  To  the  authors’ 
knowledge  only  Jin  and  Xue  [20]  have  presented  a  validated 
numerical  tool  for  the  simulation  of  SORFCs.  (Ni  et  al.  [21,22] 
address  the  modelling  and  validation  of  a  reversible  solid  oxide 
fuel  cell,  but  in  fact  they  only  report  the  SOEC  behaviour.)  If 
properly  formulated,  a  CFD  model  for  either  SOFCs  or  SOECs 
should  produce  suitable  results  in  the  other  regime  without  any 
modification  to  the  laws  (submodels)  for  the  underlying  funda¬ 
mental  physics.  However,  this  is  not  the  case  for  most  models,  as 
shown  below.  This  paper  reviews  the  challenges  posed  by  the 
development  of  such  a  unified  model,  and  how  to  overcome  them. 

Any  comprehensive  SOFC,  SOEC  or  SORFC  model  relies  to  a 
certain  extent  on  data  fitting  to  find  values  for  some  of  the  physical 
parameters,  as  it  is  the  case  of  the  charge  transfer  coefficients 
(«b,a,«fa,ab,c,af,c)  and  the  exchange  current  densities  (i0,a,i0,c)  in 
the  electrochemical  model.  This  fitting,  when  properly  resorted  to, 
is  the  consequence  of  the  incomplete  knowledge  of  the  underlying 
physical  phenomena,  or  of  the  excessive  complexity  of  such 
phenomena  for  them  to  be  accommodated  within  a  fluid-flow 
model.  The  disparity  of  spatial  scales  between  these  phenomena 
and  the  device  to  be  simulated  is  often  one  of  the  sources  of  this 
complexity. 

However,  often  this  fitting  process  is  only  vaguely  discussed 
in  the  literature.  Table  1  summarises  the  common  practices  in 
the  SOFC,  SOEC  and  SORFC  modelling  literature  to  evaluate  the 


Table  1 

Review  of  the  evaluation  of  the  electrochemical  parameters  in  SOFC,  SOEC  and  SORFC  modelling  literature./(o)  means  “calculated  as  a  function  of  o". 


Model 

Fitting 

ab,a  /  af,a 

ab,c/af,c 

io,a  /  lo,c 

Value  source 

SOFC 

[25] 

No 

0. 5/0.5 

0.5/0.5 

f(T) 

From  [56]/  guessed 

[57] 

No 

0. 5/0.5 

0.5/0.5 

5300/2000  A  m  2 

From  [58] 

[50] 

No 

0. 5/0.5 

0.5/0.5 

5300/2300  Am-2 

Not  justified 

[12] 

No 

0. 5/0.5 

0.5/0.5 

7460/10090  Am-2 

Not  justified 

[59] 

No 

0. 5/0.5 

0.5/0.5 

5300/2000  A  m-2 

From  [57] 

[60,15] 

Yes 

0. 5/0.5 

0.5/0.5 

Both  f(T,  species) 

Fitting 

[16] 

Yes 

0. 5/0.5 

0.5/0.5 

Both  f(T,  species) 

Fitting 

[45] 

Yes 

1. 5/0.5 

0.75/0.25 

Both  f(T,  species) 

Fitting 

[28] 

Yes 

2/1  [23] 

1. 5/0.5  [24] 

Both  f(T,  species) 

Fitting 

SOEC 

[61,62] 

No 

0. 5/0.5 

0.5/0.5 

Both  from  [25] 

SOFC 

[29] 

Yes 

0. 5/0.5 

0.5/0.5 

Kinetics/fitted 

Poor  fitting 

[63] 

No 

- 

- 

- 

Not  mentioned 

[21,17,26] 

No 

0.5/0.5 

0.5/0.5 

2000/5300  Am-2 

From  [57] 

[22] 

Yes 

0. 5/0.5 

0.5/0.5 

200  Am-2/ 

Not  mentioned 

[19] 

Yes 

0. 5/0.5 

0.5/0.5 

Fitted 

Poor  agreement 

[51] 

No 

0. 5/0.5 

0.5/0.5 

Both  f(T,  species) 

Literature 

SORFC 

[20] 

Yes 

0.5/0.5 

0.5/0.5 

1/0.1  Am-2 

Fitting 
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Nomenclature 

a  exponent,  Eq.  (52),  dimensionless 

A  constant  in  the  ionic  conductivity  equation,  Eq.  (A.19), 

AV-1  m-1 

ala  JANAF  constant  for  species  a,  dimensionless 

a2a  JANAF  constant  for  species  a,  I<-1 

a3a  JANAF  constant  for  species  a,  I<-2 

a4„  JANAF  constant  for  species  a,  K  3 

a5„  JANAF  constant  for  species  a,  l<  4 

b  exponent,  Eq.  (52),  dimensionless 

B  constant  in  the  ionic  conductivity  equation, 

Eq.  (A.19),  K 

B0  porous-medium  permeability-coefficient,  m2 

c  exponent,  Eq.  (53),  dimensionless 

ca  molar  density  of  species  a,  kmol  m-3 

Cp  specific  heat  at  constant  pressure  of  the  fluid, 

m2  s-2  K  1 

Cpa  specific  heat  at  constant  pressure  of  species  a, 
m2s-2I<-1 

dp  mean  pore  diameter,  m 

Va/j  ordinary  diffusion  coefficient  of  species  a  in  species  [), 

2  —  1 

m  s 

Dam  ordinary  diffusion  coefficient  of  species  a  in  a 
multicomponent  mixture,  m2  s-1 
VKa  Knudsen-diffusion  coefficient  of  species  a ,  m2  s-1 
Fact  activation  energy,  kg  m2  s-2  kmol-1 
Eb  black-body  emissive  power,  kg  s-3 

£  electric  potential  difference,  V 

£°  standard  electric-potential  difference,  V 

F  Faraday's  constant,  A  s  kmol-1 

Fj_j  view  factor  between  the  i-th  and  the  j-th  surface 
elements,  dimensionless 
h  sensible  enthalpy  of  the  fluid,  m2  s-2 

h„  sensible  enthalpy  of  species  a,  m2  s-2 

Ha  external  irradiation,  kgs-3 

i  current  density  magnitude,  A  m-2 

i  current  density,  Am-2 

i0  exchange  current  density,  Am-2 

j  a  mass  diffusion-flux  of  species  a  relative  to  the 
mass-average  velocity,  kg  m2  s-1 
( j  a)  superficial  mass  diffusion-flux  of  species  a  through  a 
porous  media,  kg  m2  s-1 
kB  Boltzmann  constant,  kg  m2  s-2  I<-1 

n  number  of  electrons  released  during  the  ionisation 

process  of  one  fuel  molecule,  dimensionless 
N  „  total  molar  flux  of  species  a,  kmol  m-2  s-1 

p  pressure,  kg  m-1  s-2 

pa  partial  pressure  of  species  a,  kg  m-1  s-2 

q  energy  flux,  kg  s-3 

Q  volumetric  heat  sources/sinks,  Joule  heat,  kg  m-1  s-3 

r  length  from  the  i-th  to  the  j-th  surface  elements,  m 

R  ideal  gas  constant,  kg  m2  s-2  kmol-1  1<-1 

Sya  mass  source  or  sink  of  species  a,  kg  m3  s_1 
Soka  Sutherland-law  parameter  for  species  a,  K 

S0fl„  Sutherland-law  parameter  for  species  a,  K 

SXo  molar  source  or  sink  of  species  «,  kmol  m3  s  -  1 
t  time,  s 

T  temperature,  K 

T  parameter  defined  in  Eq.  (A.7),  dimensionless 

T0ia  Sutherland-law  parameter  for  species  a,  K 

T opa  Sutherland-law  parameter  for  species  a,  K 

v  fluid  mass-averaged  velocity,  ms-1 


(v)  superficial  permeation  velocity,  m  s_1 
V  voltage,  V 

W  molecular  weight  of  the  fluid,  kg  kmol-1 

VV„  molecular  weight  of  species  a ,  kg  kmol - 1 

xa  molar  fraction  of  species  a,  dimensionless 

ya  mass  fraction  of  species  a,  dimensionless 

Greek  letters 

a  species,  dimensionless 

ab  backward  transfer  coefficient,  dimensionless 

af  forward  transfer  coefficient,  dimensionless 

y  pre-exponential  coefficient,  A  m-  2 

ra  dusty-gas  model  parameter,  Eq.  (12),  m2  s_1 

A*  dusty-gas  model  parameter,  Eq.  (12),  kg-1  s1  kmol 

A Hr  specific  enthalpy  of  reaction,  kg  m2  s~2  kmol-1 

e  porosity  of  the  porous  medium,  dimensionless 

ea  characteristic  energy  of  species  a  ,  kg  m2  s-2 

r/act  activation  overpotential,  V 

'Icon  concentration  overpotential,  V 

'/ohm  Ohmic  overpotential,  V 

0  angle, 

X  thermal  conductivity  of  the  fluid,  kg  m  s-3  K-1 

Xa  thermal  conductivity  of  species  a,  kgms-3  K-1 

X0a  Sutherland-law  parameter  for  species  a,  kg  ms-3  K-1 
p  fluid  viscosity,  kgm-1  s_1 

pa  viscosity  of  species  a,  kg  m-1  s_1 

p0a  Sutherland-law  parameter  for  species  a,  kgm-1  s_1 

v  number  of  occurrences  of  the  rate  determining  step  in 

an  overall  reaction,  dimensionless 
£  surface  emissivity,  dimensionless 

p  mass  density  of  the  fluid,  kg  m~3 

p  charge  density,  Cm-3 

a  ionic  conductivity,  A  V-1  m-1 

a  a  characteristic  length  of  species  a,  A 

oap  average  collision  diameter  of  species  a  and  /l,  A 
(jsb  Stefan-Boltzmann  constant,  kg  s-3  K-4 

z  tortuosity  factor  of  the  porous  medium,  dimensionless 

~vp  dusty-gas  model  parameter,  Eq.  (13),  m  s~3 

va  dusty-gas  model  parameter,  Eq.  (13), 

kg-1  m-1  s1  kmol 

if  dusty-gas  model  parameter,  Eq.  (14),  m  s_1 

~&a  dusty-gas  model  parameter,  Eq.  (14), 

kg-1  m-1  s1  kmol 
cj)  electric  potential,  V 

</>a  electric  potential  of  the  electric-conducting  phase  of 

the  fuel  electrode,  V 

<j>c  electric  potential  of  the  electric-conducting  phase  of 

the  air  electrode,  V 

<pe  electric  potential  of  the  ion-conducting  phase  of  the 

electrolyte,  V 

<p*  effective  electric  potential  of  the  ion-conducting  phase 

of  the  electrolyte,  V 

@a/i  semi-empirical  correlation  parameter, 

Eq.  (A.12),  dimensionless 
a >  heat  of  reaction,  kg  m-1  s-3 

£2d  collision  integral,  dimensionless 

Subscripts 

a  fuel  electrode 

b  bulk,  feeding  mixture  concentration 


704 


M.  Garcia-CampruM  et  al.  /  Renewable  and  Sustainable  Energy  Reviews  33  (2014)  701-718 


c 

air  electrode 

rc 

cathode  reaction  wall,  i.e.  cathode-electrolyte 

e 

electrolyte 

interface 

eq 

equilibrium 

rad 

radiation 

f 

fluid 

ref 

reference  value 

i 

isothermal-surface  index 

rev 

reversible 

ia 

interconnect/fuel-electrode  interface 

s 

solid 

ic 

interconnect/air-electrode  interface 

SOEC 

at  SOEC  mode 

irrev 

j 

irreversible 

isothermal-surface  index 

SOFC 

at  SOFC  mode 

oc 

r 

open  circuit 

generic  reaction  wall,  i.e.  electrode-electrolyte 

Superscripts 

ra 

interface 

anode  reaction  wall,  i.e.  anode-electrolyte  interface 

eff 

effective  value  in  the  porous  media 

electrochemical  parameters.  It  is  widely  assumed  that  the  charge 
transfer  coefficients  are  equal  to  0.5,  what  would  imply  a  single  step, 
single  electron  reaction.  Theoretically,  it  is  clear  that  this  assumption 
does  not  apply  to  the  hydrogen  reduction/oxidation  reaction.  Experi¬ 
mentally,  there  are  interesting  works  refuting  this  hypothesis  [23,24], 
Numerically,  there  is  also  an  interesting  contribution  of  Grondin  et  al. 

[19]  stating  that  Butler-Volmer  equation  together  with  the  single 
step,  single  electron  hypothesis  does  not  permit  relevant  computing 
predictions  of  SOEC  polarisation  curves.  However,  this  assumption  is 
still  used  as  it  is  simple  and  provides  (apparently)  good  fitting  to 
experimental  data.  In  this  work,  we  show  that  this  assumption  fails  if 
the  model  is  properly  validated  and  that  more  accurate  values  can  be 
found  by  means  of  a  thorough,  careful  fitting  process. 

Once  the  electrochemical  reaction  is  assumed  as  a  single  step, 
single  electron  one,  the  exchange  current  densities  are  the 
remaining  unknowns  of  the  model.  How  these  are  found  and 
reported  in  the  published  literature  deserves  some  close  scrutiny. 
As  shown  in  Table  1,  some  papers  in  the  literature  do  not  mention 
these  values.  Sometimes,  when  no  experimental  data  is  found 
these  values  are  guessed  without  any  further  validation  [25],  Some 
authors  take  constant  values  from  published  experimental  works, 
even  if  the  operating  conditions  or  materials  of  the  experimental 
cell  do  not  match  those  of  the  cell  being  modelled  [21,26],  For 
instance,  Ni  et  al.  [26]  assume  that  the  exchange  current  density  of 
the  carbon  dioxide  electrolysis  is  that  of  the  water  electrolysis.  Ni 
et  al.  [21]  validate  the  numerical  data  with  the  experimental  data 
of  Momma  et  al.  [27],  at  first  sight  this  validation  would  indicate 
that  taking  constant  values  from  the  literature  may  lead  to 
reasonable  fitting  with  experimental  data.  However,  close  exam¬ 
ination  of  the  experimental  data  used  [27]  reveals  that  only  the 
linear  range  of  the  experimental  characteristic  curve  is  used 
(ie(0, 5000)  Am-2),  avoiding  the  comparison  with  the  experi¬ 
mental  data  at  high  current  densities  where  the  experimental 
trend  is  not  linear  anymore.  Table  1  also  shows  that  the  value  of 
the  exchange  current  densities  is  estimated,  in  other  works,  to  fit 
the  numerical  results  to  an  experimental  data  set.  For  example,  in 

[20]  a  constant  value  for  both  exchange  current  densities  is 
guessed  to  fit  the  experimental  data,  neglecting  its  dependence 
on  temperature  and  species  concentration;  the  resulting  values  are 
unusually  small  (1,0.1  Am-2).  This  is  also  done  in  [22]  for  one 
exchange  current  density,  but  the  authors  do  not  indicate  whether 
it  is  the  cathodic  or  the  anodic  one.  In  the  SOFC  literature,  it  is  a 
common  practice  to  use  experimental  correlations  that  express  the 
exchange  current  densities  as  a  function  of  the  temperature  and 
species  distribution  at  the  electrodes,  e.g.  [28],  These  are  generally 
Arrhenius  type  expressions,  the  pre-exponential  factors  of  which 
are  the  fitted  parameters.  In  the  SOEC  literature  this  kind  of 
correlation  is  used  in  [29]  leading  to  a  very  poor  agreement. 

From  the  above  review,  the  authors  conclude  that  the  evalua¬ 
tion  of  the  electrochemical  parameters  is,  in  general,  not  properly 


justified.  In  our  opinion,  this  not  only  greatly  diminishes  the 
usefulness  of  modelling  as  a  design  tool,  but  also  contributes  with 
uncertainty  rather  than  guidance  when  published  in  the  literature. 

In  the  present  work  we  introduce  a  comprehensive  model  for 
the  simulation  of  an  SORFC,  together  with  a  detailed  description  of 
how  the  electrochemical  parameters  are  fitted.  The  validity  of  the 
proposed  model  and  fitted  parameters  is  also  shown  by  compar¬ 
ison  of  the  numerical  results  with  experimental  data  for  a  micro¬ 
tubular  SORFC  [10], 

The  model  presented  caters  for  the  following  phenomena: 
(a)  mass  transport  through  channels;  (b)  momentum  transport 
through  channels  and  porous  solids;  (c)  multispecies  mass  transfer 
through  porous  media  (convection,  ordinary  diffusion,  Knudsen 
diffusion);  (d)  heat  transfer  in  gases,  porous  media  and  solids 
by  conduction,  convection  and  radiation;  (e)  charge  transport 
through  an  impervious  solid;  and  (f)  reversible  electrochemistry. 
The  model  has  been  implemented  in  OpenFOAM-1.5-dev  [30],  an 
open-source  CFD  platform  based  on  the  finite-volume  method. 
The  full  mathematical  model  and  the  relevant  features  of  the 
numerical  method  are  thoroughly  described  in  Sections  2  and  3 
respectively.  The  model  can  be  regarded  as  a  compendium  of  the 
most  appropriate  submodels  found  in  the  literature  for  the  task; 
and  the  full  description  of  this  comprehensive  model  in  this  paper 
should  enable  its  implementation  in  any  code  by  any  interested 
researcher. 

The  model  is  validated  using  experimental  data  presented  in 
Section  4.  The  unavoidable  fitting  process  is  next  described  in 
Section  5.  It  reveals  that  some  common  simplifications  used  in 
the  SOFC  and  SOEC  modelling  literature  (see  Table  1)  might  be 
inappropriate.  Our  findings  highlight  the  often-ignored  relevance 
of  this  fitting  process  for  obtaining  a  predictive  tool  that  can  be 
used  reliably  across  the  range  of  SOFC,  SOEC  and  SORFC  designs. 

Finally  in  Section  6  our  numerical  results  are  compared  with 
experimental  data  to  validate  the  model.  Not  only  the  numerical 
and  experimental  characteristic  i-V  curves  agree  well,  but  so  do 
the  magnitudes  of  other  variables  that  depend  on  the  fitted 
parameters,  such  as  the  exchange  current  density  and  the  activa¬ 
tion  overpotential. 


2.  Mathematical  model 

In  this  section  the  mathematical  model  that  describes  the 
operating  principles  of  an  SORFC  is  presented.  The  model  is  an 
extension  of  one  developed  previously  for  SOFC  simulation  [16,31], 
the  new  features  being:  (i)  a  new  multidimensional  charge- 
transfer  model  for  the  electrolyte;  (ii)  an  improved  formulation 
of  the  calculation  of  heat  sources  and  sinks;  (iii)  the  extension  of 
the  electrochemical  model  to  SOEC  conditions. 
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The  operating  principle  of  a  solid  oxide  regenerative  fuel  cell,  both 
in  SOFC  and  SOEC  modes,  involves  a  complex  network  of  interde¬ 
pendent  physical  phenomena  within  the  several  SORFC  components. 
The  governing  equations  that  describe  a  given  physical  phenomenon 
may  change  from  one  component  to  another  because  of  the 
components’  differing  nature,  e.g.  porous  media,  gas  channels  or 
impervious  solids.  To  describe  each  phenomenon  with  the  appro¬ 
priate  equations,  the  mathematical  model  has  been  split  into  a  set  of 
submodels,  namely:  (i)  channel;  (ii)  electrode;  and  (iii)  electrolyte 
submodels.  Two  additional  submodels  account  for  the  radiative  heat 
transfer  and  the  electrochemistry;  these  operate  at  the  interface 
between  some  of  the  components  (respectively  channel-electrode 
and  electrode-electrolyte). 


2.1.  Channel  model 


The  SORFC  channels  are  ducts  along  which  the  reactant 
mixtures  are  driven  from  the  SORFC  inlet  to  the  electrodes.  The 
byproducts  of  the  electrochemical  reactions  exit  the  electrode  and 
join  the  unconverted  reactants  in  the  channel,  where  they  are 
evacuated  to  the  SORFC  outlet.  The  channel  model  accounts  for  the 
phenomena  summarised  in  Table  2  and  described  below. 

The  differential  form  of  the  continuity  equation  is 


^+V-0f)  =  0  (1) 

ot 

where  v  is  the  fluid  velocity  and  p  is  the  fluid  density,  given  by 
the  ideal  gas  law; 


P  = 


pW 
~ RT 


(2) 


where  p  is  the  pressure,  W  the  mixture  molecular  weight,  T  the 
temperature  and  R  the  ideal  gas  constant. 

The  momentum-conservation  equation  for  a  Newtonian  fluid 
with  negligible  bulk  viscosity  (pv  ~  0),  with  negligible  body  forces, 
is  (in  vector  form) 

dip  V  )  — >• — >-  — >. 

-'C-2+V  (/jV  v)-V-fl(VV)  =  V- 
at 

(3) 


MV'v)T-^tr[(Vv’)T]  I 


-Vp 


We  model  j  a  according  to  the  consistent  effective  binary 
diffusion  method  [32]: 

j  a  —  —  pC>amSy a+  py aYSDfimSy fi)  (5) 

V/l 


where  Vam  is  the  diffusion  coefficient  of  species  a  in  a  multi- 
component  gas  mixture  (see  Appendix  A.2).  This  formulation 
ensures  that  the  diffusive  fluxes  relative  to  averaged  fluid  velocity 
properly  add  up  to  zero.  It  thus  combines  the  simplicity  of  the 
effective  binary  diffusion  approach  with  flux  consistency. 

From  Eqs.  (4)  and  (5),  the  conservation  equation  of  a  chemical 
species  a  in  the  SORFC  channel  is 


d(pya) 

dt 


+  V  •  (pya~v  )- V  •  (pVamVya)  + V  • 


pya'DPpmNyp 

vp 


=  s, 


(6) 


The  conservation  of  energy  in  the  system  is  accounted  for 
through  the  following  equation  for  the  sensible  enthalpy: 

pCp^+  V  ■  (pCp~v  T)  -  V  ■  (AVT)  =  -  I(yA)^+  TV 

OL  ya  OL 


■  (.Cpp~v  )  —  V  • 


K  j  A) 

.Vo 


(7) 


where  Cp  is  the  fluid  specific  heat  at  constant  pressure 
(see  Appendix  A.5),  A  is  the  thermal  conductivity  of  the  fluid 
(see  Appendix  A.4),  and  ha  is  the  sensible  enthalpy  of  species  a 
(see  Appendix  A.6). 

In  Eq.  (7)  the  following  assumptions  and  provisions  have  been 
made  [16]:  (i)  the  sensible  enthalpy  is  defined  as  dh  =  CpdT ; 
(ii)  the  specific  heat  at  constant  pressure  depends  on  temperature 
Cp=f(T)  (see  Appendix  A.5);  (iii)  the  diffusion  of  heat  is  due  to 
both  heat  conduction  (Fourier’s  law)  and  species  mass  diffusion; 

(iv)  the  power  delivered  to  the  fluid  by  body  forces  is  negligible; 

(v)  viscous  dissipation  of  energy  can  be  neglected;  (vi)  large 
pressure  differences  in  the  system  are  not  expected;  (vii)  the 
reaction  heat-release  within  the  channel  is  null  since  the  reaction 
sites  are  located  in  the  electrodes;  and  (viii)  radiation  is  from 
surface  to  surface,  and  hence  there  are  no  volumetric  sources  of 
radiative  heat  in  the  channels  [33,34], 


2.2.  The  electrode  model 


where  p  is  the  dynamic  viscosity  of  the  fluid  (see  Appendix  A.l). 

The  differential  form  of  the  equation  for  the  conservation  of 
chemical  species  a  in  the  channel  is 

^+V.^)  +  V.r„=S,  (4) 

where  ya  is  the  mass  fraction  of  species  a,  j  a  is  the  mass 
diffusion-flux  of  species  a  relative  to  the  mass-average  velocity, 
and  Sya  stands  for  the  volumetric  sources  or  sinks. 


Table  2 

Summary  of  the  main  features  of  the  mathematical  submodels. 


Submodel 

Phenomena 

State/constitutive  equations 

Channel 

Mass  transport 

Ideal  gas 

Momentum  transport 

Newtonian  fluid 

Species  transport 

Consistent  diffusion 

Energy  transport 

Fourier  multicomponent 

Electrode 

Mass  transport 

Ideal  gas 

Momentum  transport 

Darcy 

Species  transport 

Dusty-gas 

Energy  transport 

Fourier  multicomponent 

Electrolyte 

Charge  transport 

Ohm 

Energy  transport 

Fourier 

Radiation 

Surface-to-surface 

View-factor  method 

Electrochemical 

Redox 

Nernst,  Butler-Volmer 

The  electrode  model  accounts  for  the  momentum,  mass,  species 
and  heat  transport  within  the  SORFC  electrodes;  both  anode  and 
cathode  are  porous  ceramics  that  are  assumed  to  be  homogeneous 
and  isotropic. 

The  conservation  equation  for  the  species  a  in  a  porous 
medium,  i.e.  the  counterpart  in  the  electrode  to  Eq.  (4)  in  the 
channel,  is 


v  •  (pya(v  » + v  •  (fa)  =  sya  (8) 

where  (  v  )  is  the  fluid  permeation  velocity  through  the  electrode 
(or  superficial  velocity)  and  (  j  a)=  e  j  a  represents  the  superficial 
combined  diffusion  flux  (ordinary  and  Knudsen)  through  the 
porous  medium. 

The  dusty  gas  model  (DGM)  is  the  most  convenient  approach  to 
model  such  a  species  transfer,  as  reported  in  [35,36],  It  is  worth 
noting,  though,  that  the  theoretical  basis  of  the  DGM  [37]  was 
criticised  by  Kerkhof  [38,39],  who  as  an  alternative  proposed  the 
binary-friction  model  (BFM).  However,  the  BFM  does  not  improve 
the  DGM  results  [38]  and  both  models  showed  similar  agreements 
and  discrepancies  with  experimental  data.  The  DGM  is  thus  used 
in  this  work;  its  general  form  is  as  follows  [37]: 


-JLvn  ^  y,  XpNa-xaNp  ,  N 

p  t  v  ” a  ^  _ eff 

Rl  T>£, 


V 


0 


+ 


1  Pa  Bo 

tfffRTp 


Vp 


(9) 
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where  N„  is  the  total  molar  flux  of  species  a,  B0  is  a  constant 
porous-medium  permeability-coefficient,  and  =  (e/z)T>al!  and 
=  ( e/T)VKa  stand  for  the  effective  binary  and  Knudsen  diffu¬ 
sion  coefficients  respectively  (see  Appendices  A.2  and  A.3). 

Considering  the  total  molar  flux  of  species  a,  Na  = 
W-\pya(v)+(Ta ».  Eq.  (8)  is  rewritten  as 


d(eC„) 

dt 


+  V  ■  N  a— S, 


(10) 


where  Sxa  is  the  volumetric  molar  sources  or  sinks  of  species  a, 
SXa  =  W;'Sya. 

Replacing  Eq.  (9)  into  Eq.  (10)  and  rearranging  it  results  in 
v  ■  (Clvpj  +  v  •  (vfpJ  +  V-  (Iff  PJ  =  SXa  (11) 


where  the  following  coefficients,  here  referred  to  as  DGM  para¬ 
meters,  have  been  defined  for  compactness  and  ease  of  interpreta¬ 
tion: 


r 


a 


l 


(12) 


(13) 


0  a 

RT 


va=r„RTZ 

P*a 


(14) 


ra  (m2  s_1)  represents  a  global  diffusion  coefficient  of  species  a  in 
the  multicomponent  gas  mixture  in  a  porous  media  ;  ~vpa  (m  s~]) 
is  a  superficial  velocity  due  to  the  pressure  gradient  in  the  porous 
medium;  and  ~o  a  (m  s_1)  is  a  superficial  velocity  induced  by  the 
flow  of  the  other  species. 

Eq.  (11),  applied  to  all  the  components  of  the  gas  mixture, 
represents  the  complete  set  of  equations  for  mass  transport,  and 
do  not  need  to  be  supplemented  by  any  additional  equations  of 
motion  [37],  The  global  pressure  p  can  be  computed  as 


P=lP, 

Va 


(15) 


The  diffusive  heat  fluxes  appearing  in  Eq.  (17)  are  given  by 


Q  f  —  ~ 2/VT  +  2(  j  aha) 

Va 


(18) 


~Q  s=  — 2SVT  (19) 

An  effective  thermal  conductivity  of  the  porous  medium  can  be 
defined  as  Aett  =  eAt  +  (1  -e)As  [43].  Considering  Eqs.  (18),  (19)  and 
the  following  relationship: 

Pfhf(v)+X((Ta)ha)  =  I/[ha(Wa~Na)]  (20) 

Va  Va 


then  Eq.  (17)  can  be  rewritten,  after  some  rearrangement,  as 
dT 


[ePfCpf  +  (1  —  e)ftCp,s]-^:+  V 
=  TV  • 


T-Z(CPaWaNa) 

.  Va 


-  V  •  (Aeff  VT) 


X(CpaWaN0 

Va 


-2  [haV-(WaNa)] 

Va 


dpf 


dps 


-fhr — E_(l  -e)hs-LA-ym 


dt 


dt 


(21) 


where  w,  the  heat  of  reaction,  has  two  contributions:  (a)  due  to  the 
reversible  reaction  heat  release  ( Qj-ev ),  this  being  released  only  in 
the  fuel  electrode;  and  (b)  due  to  the  irreversibilities  of  the  process 
(Qin-ev),  this  being  released  where  such  irreversibilities  or  over¬ 
potentials  take  place.  Thus, 

f  Qirrev.c  in  the  air  electrode 

m  ~  1  Q_rev  +  Qirrev.a  in  the  fuel  electrode 


The  chemical  reaction  is  assumed  to  take  place  on  the  electrode¬ 
electrolyte  interfaces,  and  thus  this  is  where  the  heat  of  reaction  is 
released;  it  is  therefore  accounted  for  in  Eq.  (21)  as  a  heat  flux 
boundary  condition  (see  Section  3). 


2.3.  The  electrolyte  model 

The  electrolyte  is  an  impervious  solid  traversed  by  an  anionic  (02~ ) 
flux.  The  model  consists  of  equations  for  the  conservation  of  charge, 
which  simultaneously  enforces  mass  conservation,  and  of  energy. 

The  charge-conservation  equation  is,  considering  Ohm's  Law,  as 
follows: 

^-V-(<feV<fe)  =  0  (23) 


and  the  superficial  permeation  velocity  may  then  be  estimated 
using  Darcy’s  Law  [40]: 

(v)=  Vp  (16) 

P 

The  authors  have  published  an  open-source  species  mass  transfer 
library  that  includes  this  DGM  algorithm,  it  is  available  in  [41  ]. 

The  sensible  enthalpy  conservation  equation  in  a  porous 
medium,  assuming  local  thermal  equilibrium  between  the  porous 
matrix  and  the  fluid  flowing  through  the  interstitial  space  [42],  is 

ahp//t/  +  n-g)pshs]  +  v  +  [t.7f/  +  (1  _e)lfs]  =  m  (17) 

where  the  subscripts /and  s  stand  for  fluid  and  solid  respectively, 
and  the  following  assumptions  and  provisions  have  been  made: 
the  power  delivered  to  the  fluid  by  the  body  forces  and  the  viscous 
dissipation  are  negligible;  there  are  no  large  pressure  differences 
in  the  system;  and  there  are  no  volumetric  heat  sources  within 
the  electrodes  because  Joule  heating  takes  place  only  within  the 
electrolyte,  and  radiation  can  be  neglected  in  the  opaque  electro¬ 
des  [34,33],  Note  that  surface-to-surface  radiative  heat  flux,  q  md, 
is  considered  at  the  electrode  outer  boundary  by  means  of  a 
boundary  condition  (see  Section  2.4). 


where  p  is  the  volumetric  charge  density,  de  is  the  ionic  con¬ 
ductivity  of  the  electrolyte  (see  Appendix  A.7),  and  tpe  is  the 
electric  potential  of  the  ion-conducting  phase  of  the  electrolyte. 

Heat  transfer  in  the  electrolyte  is  due  to  conduction.  Considering: 
(i)  Fourier’s  law,  q  =  —XS/T\  (ii)  Joule  heating,  Q  =  f2  cfe  1;  and 
(iii)  the  definition  of  sensible  enthalpy,  dh  =  Cp  dT,  then  the  energy 
conservation  equation  can  be  written  in  terms  of  temperature  as 

pCp V  •  (AVT)  =  -hCK'  (24) 

dt  dt  ae 

where  i  =  \  i  |  is  the  current  density  magnitude. 

For  those  cell  configurations  where  the  electrolyte  is  exposed 
to  the  channels,  such  as  the  case  of  a  microtubular  cell  where 
the  cathode  is  shorter  than  the  electrolyte,  the  surface-to-surface 
radiative  heat  exchange  is  included  as  a  heat-flux  boundary  condition 
using  a  surface-to-surface  radiation-model  (see  Section  2.4). 

2.4.  The  radiation  model 

The  importance  of  radiative  heat  transfer  has  been  discussed  in 
the  SOFC-modelling  literature  [33,34,14].  The  prevailing  ideas  are: 
(i)  the  most  common  gases  in  SOFCs  behave  as  non-participating 
media,  and  thus  surface-to-surface  radiation  exchange  is  the 
only  radiative  transfer  mode  that  needs  to  be  considered  in  the 
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channels  of  an  SOFC;  (ii)  the  electrodes  are  opaque  to  radiation; 
(iii)  the  effect  of  volumetric  radiation  within  a  thin  electrolyte  is 
negligible,  though  the  electrolyte  material  is  optically  thin.  These 
conclusions  are  also  applicable  to  SOEC. 

The  accurate  calculation  of  surface-to-surface  radiation  requires  of 
complex  models  accounting  for  both  parallel  and  oblique  radiation 
[14],  such  the  view-factor  method  used  in  the  present  work  [44]: 


Qrad.i  =  ^i  \  Eb.i  ~  Z(f  i -jEb j)  ~  Ho,i  +  Z 

i  i 


(25) 


where  qrad  i  is  the  surface-to-surface  radiative  heat  flux  arriving  at  the 
infinitesimal,  isothermal  surface  element  dA,;  and  j  indicates  an 
infinitesimal,  isothermal  element  of  surface  dAj  irradiating  element  i. 
In  Eq.  (25),  f  is  the  emissivity  of  the  surface,  Eb  =  o^T4  is  the 
blackbody  emissive  power,  asb  being  the  Stefan-Boltzmann  constant; 
H0  is  the  incident  radiation  entering  or  leaving  the  domain  through  an 
opening  (H0= 0,  in  this  work),  and  F,-_ j  is  the  view  factor  between  the 
two  infinitesimal  elements  dA,  and  dAj.  The  view  factors  are  defined  as 


cos  d,  cos  0:  ,, 
Fi-i  = - ^ - J-dA< 


(26) 


where  d,  and  d,  are  the  angles  between  the  normal  vector  to  the 
surface  elements  dA,  and  dA,  and  the  line  (of  length  r)  connecting 
them.  Further  details  of  the  radiation  model  are  given  elsewhere  [16], 


2.5.  The  electrochemical  model 


The  electrochemical  model  represents  the  electrochemical 
reaction  that  takes  place  at  the  triple  phase  boundaries  (TPBs), 
where  electrolyte,  reactants  and  the  catalyst  region  of  the  elec¬ 
trode  meet.  These  TPBs  are  confined  to  a  region  close  to  the 
electrode-electrolyte  interface  and  therefore  the  reaction  is  repre¬ 
sented  as  occurring  on  a  surface. 

In  the  following  subsections  the  electrochemical  models  for  the 
SOFC  and  SOEC  operating  modes  are  presented. 


2.5.1.  SOFC  mode 

The  electrochemical  reaction  that  takes  place  in  a  hydrogen-fed 
SOFC  is 

h2+|o2-*h2o 

The  standard  expression  for  the  cell  voltage  of  such  a  reacting 
system  is 

^  ==  Voc  —  dconc.a  ~  dconc.c  dohm  ~  1 1  act. a  ~  dact.c  (27) 

where  subscripts  a  and  c  stand  for  the  fuel  and  air  electrodes 
respectively  (see  Table  3)  and  yconc,  r/ohm,  r/act  are  the  concentration, 
ohmic  and  activation  overpotentials  respectively.  These  overpo¬ 
tentials,  which  are  positive,  are  defined  below: 


Concentration  overpotential: 


RT. 

d cone,  a  —  ' 


PH2.ttPH2O.ra 

PH2O.bPH2.ra 


(28) 


RT , 

dconc.c  —  2 ~p  ^ 


(29) 


Table  3 

Electrochemical  roles  of  the  cell  electrodes  for  both  operating  modes. 


Electrode 

Label 

SOFC  mode 

SOE  mode 

Fuel  electrode 

a 

Oxidation 

Reduction 

Air  electrode 

c 

Reduction 

Oxidation 

•  Ohmic  overpotential:  From  Eq.  (27) 

Wohm  =  V oc  —V  —  iJCOnc,a  ~  Vconct  ~  Vact^a  ~  ^actx 

where  by  definition  [45] 

*7act,a  =  Ea  ~  Eeq,a 

Vactx  =  Eeq,c  ~  Ec 


Voc  —  V rev  ~E  *1  cone, a  ~E  *7conc,c 

with  Eeqa,  Eeqx ,  Ea,  Ec  and  Vrev  being  defined  as 
Ph2,m 


Eeq,a  — 


Pn2o,r 


E 


eq,c  — 


RT 
2 F 


In 


1/21 


(30) 

(31) 

(32) 

(33) 

(34) 

(35) 


Ea  =  4>a  ~  4>e.a  (36) 

Ec  =  <Pc~  4>e.c  (37) 


V  rev  —  Eeq.c  —  Eeq.a  (38) 

Here  rpa  and  tpc  are  the  electric  potentials  of  the  electron¬ 
conducting  phase  in  the  fuel  and  air  electrodes;  and  rpea  and 
4>e.c  are  the  electric  potentials  of  the  ion-conducting  phase  in 
the  fuel  and  air  electrodes  respectively. 

From  Eqs.  (31)  to  (38),  Eq.  (30)  may  be  rewritten  as 

dohm  =  tpe.ra  ~  fie.rc  T  fic.rc  ~  (l'a.ca  "  ^  (39) 

In  the  present  work,  the  electrodes  are  assumed  to  be  ideal 
electron  conductors: 

Vi/)a  =  0  =>  a  =  constant  =>  <pam  =  tpaia  (40) 

Vr/>c  =  0  =>•  rj>c  =  constant  =>•  tpcrc  =  ipcic  (41) 

where  <paia  and  <pcic  are  the  electric  potential  of  the  electron- 
conductor  phase  of  the  electrode  at  the  electrode-interconnect 
interface,  that  satisfy  the  following  equation: 

17  =  4>c,ic  -  4>a.ia  (42) 

Therefore,  from  Eqs.  (40),  (41),  (42)  and  (39),  the  ohmic  over¬ 
voltage  is  defined  as 

dohm  =  tPe.ra—  fie.rc  (43) 

The  ohmic  overpotential  is,  in  Eq.  (27),  defined  as  a  positive 

magnitude;  therefore  from  Eq.  (43): 

dohm  ->  0  =$■  (pe  ra  ^  tpe.rc  (44) 


Eq.  (44)  conflicts  with  Ohm's  law  (<pe  m  <  <pe  rc )  as  the  current 
flows  from  the  cathode  to  the  anode.  To  remedy  this,  a  new 
variable  is  used: 


<P*e  =  ~<t>e 

(45) 

Eqs.  (31)  and  (32)  may  be  rewritten  as 

Uactfi  =  fiaja  ~E  tPeja  ~  Eeq,a 

(46) 

tlactx  =  Eeq,c  ~  fieje  ~  $e,rc 

(47) 

From  Eqs.  (30),  (46)  and  (47) 

tlohm  =  fiejc  ~  fieja  ~E  *Pc,rc  —  $ a,ra  ~  ^ 

(48) 

v 


Therefore,  the  ohmic  overpotential  is  redefined  as 


dohm  —  fie.rc  i’e.ra 


(49) 
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where  r/0hm  is  positive  as  required  and  4>*  satisfies  Ohm's  law 

(<ra<<rc)- 

•  Activation  overpotential: 

This  is  given  by  the  most  general  form  of  Butler-Volmer 
equation  [46]: 


ha  —  lo.a 


exp 


«b  ,aF>lac 
RT 


-exp 


l/f.aRllact..a\ 

RT  ) 


(50) 


exp 


(  abxFilactA 

V  RT  ) 


-exp 


RT  ) 


(51) 


where  ab  and  af  are  the  backward  and  forward  transfer 
coefficients,  and  i0  0  and  i0-c  are  the  exchange  current  densities 
at  each  electrode,  for  the  calculation  of  which  the  following 
experimental  correlations  are  widely  used  [47,48]: 

io,a  =  Ya  XH2  XH20  exP  ^  (^2) 


io,c  =  Yc  xo2  exP  (  - (53) 

with  Eact,a  and  £act,c  being  the  anodic  and  cathodic  activation 
energies:  a,  b  and  c  expressing  the  concentration  dependency; 
and  ya  and  yc  being  the  anodic  and  cathodic  pre-exponential 
coefficients  respectively. 


2.5.2.  SOEC  mode 

If  the  cell  operates  as  an  SOEC,  the  electrochemical  reaction  is 
reversed: 


H20^H2+i02 

The  cell  voltage  is  given  by  Eq.  (27),  where  the  open  circuit  voltage 
of  the  reverse  reaction  is  to  be  considered:  VocSOec  =  -Vocsofc- 
Hence, 


^  —  ^oc.SOEC  Pconc.a  tfconc.c  *lact,a  1 lact.c  Pohm  S  0  (54) 

The  overpotentials  are  positive  magnitudes,  and  their  definitions 
slightly  differ  from  those  in  the  SOFC  mode,  as  indicated  below: 


•  Concentration  overpotential: 


*1  cone, a  z 


—  ^Zqn  |  Rh2 o,bPhi2  ,ra 


2  F  \PH2.bPH20.ra 


(55) 


to  the  air  electrode  and  thus 

Pohm  >  0  =y  fpe.ra  ''  (Pejc  4>e  =  *Pe  (51 ) 


In  summary,  the  electrochemical  model  consists  of  a  set  of 
algebraic  equations,  namely  Eqs.  (27),  (28),  (29),  (49),  (50)  and  (51) 
for  SOFC  mode  and  Eqs.  (27),  (55),  (56),  (49),  (50)  and  (51)  for 
SOEC  model.  However,  the  transfer  coefficients  (ab,af)  and  the 
exchange  current  densities  (i0,a,io.c)  in  Eqs.  (50)  and  (51)  are  cell 
dependent  and  not  measured  experimentally.  Hence,  for  comple¬ 
tely  specifying  a  model,  these  values  are  often  guessed  in  a  data 
fitting  process;  in  this  work,  this  is  presented  in  Section  5. 


3.  Numerical  details 

The  numerical  solution  of  the  mathematical  model  introduced 
in  Section  2  has  been  addressed  by  developing  an  in-house 
algorithm  solvable  in  OpenFOAM-1.5-dev  [30]  using  a  finite- 
volume  method.  In  this  work  we  study  the  steady-state  SORFC 
operation,  and  therefore  the  time  derivatives  in  the  transport 
equations  are  null. 

To  facilitate  the  implementation  of  the  mathematical  model, 
the  spatial  domain  is  split  into  subdomains  for  each  of  the  five  cell 
components;  each  subdomain  is  meshed  separately  resulting  in 
five  conformal  adjacent  meshes  (Fig.  1),  namely:  (i)  fuelchan- 
nelMesh  for  the  fuel  channel  domain;  (ii)  anodeMesh  for  the  fuel 
electrode;  (iii)  electrolyteMesh  is  the  mesh  for  the  electrolyte 
subdomain;  (iv)  cathodeMesh  for  the  air  electrode;  and  (v) 
airchannelMesh  for  the  air  channel.  These  numerical  grids  are 
adjacent  to  each  other;  between  two  adjacent  meshes  there  is  only 
one  physical  boundary  but  this  is  represented  twice  in  the 
numerical  model,  one  for  each  of  the  two  adjacent  meshes  (these 
are  the  so-called  coupled  boundaries). 

The  overall  algorithm  linking  the  several  solution  procedures  is 
embodied  in  a  main  code.  Its  salient  algorithmic  steps  are  shown 
in  Fig.  2.  First,  the  computational  meshes  are  created  and  checked, 
returning  a  fatal  error  if  they  are  not  face-matching.  Then,  the 
fields  required  in  each  spatial  domain  are  created  and  their  values 
initialised.  Next,  the  radiation  model  is  solved  providing  the 
radiative  heat  fluxes  to  be  introduced  as  a  source  term  in  the 
energy-conservation  equation  for  the  air  electrode  and  the  elec¬ 
trolyte  [16],  Hence,  from  Eq.  (24)  the  energy-conservation  equa¬ 
tion  in  the  electrolyte,  considering  surface-to-surface  radiation 
and  steady  state,  is 


W,c-2Fln(^§)  (56) 

•  Activation  overpotential:  ijac[a  and  qactc  are  given  by  Eqs.  (50) 
and  (51)  respectively,  where  ioa  and  ioc  are  given  by  Eqs.  (52) 
and  (53),  since  the  exchange  current  density  is  by  definition 
the  same  for  the  direct  (SOFC)  and  reverse  (SOEC)  reactions; 
however  the  charge  transfer  coefficients  change  according  to 


ab,a.S0EC  =  af.a,SOFC  (57) 

ab,cSOEC  =  “f .C.SOFC  (58) 

af,a,SOEC  =  ab.a,SOFC  (59) 

af,c,SO£C  =  ab,c,SOFC  (60) 


•  Ohmic  overpotential:  riohm  is  defined  in  Eq.  (43)  as  ijohm  = 
‘t’eja  -  tlle,rc-  In  an  SOEC,  the  current  flows  from  the  fuel  electrode 


-V-(2VT)  =  ^— V-  qrad  (62) 

where  ~q  rad  is  non-zero  only  at  the  radiative  boundaries  of  the 
electrolyte.  Similarly,  from  Eq.  (21)  the  energy-conservation  equa¬ 
tion  in  the  air  electrode,  considering  surface-to-surface  radiation 
and  steady-state  conditions,  is 


V  • 


T'Z(CpaWaNa) 

.  Va 


V  •  VT) 


=  rv- 


Z(CpaWaNa ) 

.Va 


2[/t0V-(W„JVa)]-V.-qrad 

Va 


(63) 


where  q  rad  is  non-zero  only  at  the  radiative  boundaries  of  the  air 
electrode. 

The  electrode  model  is  solved  at  both  electrodes  and  the 
relevant  information,  such  as  velocity,  pressure  and  species  con¬ 
centrations  at  the  coupled  boundaries,  is  transferred  to  the 
adjacent  channel  sub-domains.  Chemical  reaction  is  modelled 
as  a  surface  reaction;  therefore,  in  Eq.  (11)  SXii  =  0  and  the  mass 
source/sink  is  imposed  at  the  reaction  wall  as  a  fixed  gradient 
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Air  channel 


Air  electrode 


Electrolyte 


Fuel  electrode 


Fuel  channel 


Fig.  1.  Schematic  of  the  numerical  multi-domain. 


SORFC-solver 

{ 

f or (startTirae ,  startTime  <  endTime,  startTime  +=deltaTime) 

{ 

create:  fuelChannelMesh,  anodeMesh,  electrolyteMesh,  cathodeMesh, 
airChannelMesh ; 
checkCoupledPatches ; 

create:  fuelChannelFields ,  anodeFields,  electrolyteFields ,  cathodeFields , 
airChannelFields ; 
f ieldslnitialization; 

radiationSolve ; 

electrodesSolve ; 
electrodeToChannelCoupling ; 


channelsSolve ; 
channelToElectrodeCoupling ; 


electrochemistrySolve ; 
electrolyteSolve ; 

checkBalances ; 
write ; 

} 

exit 

} 


boundary  condition  given  by 


Fig.  2.  Structure  of  the  segregated  solver  for  the  mathematical  model  of  an  SORFC. 

from  Eq.  (63): 

-  V  •  (Aeff  VT) 


{ 2 

Xftj  N  a,r  —  Xaj  N  jjj 

v  • 

TZ(CpawaNa) 

veff 

L  UapJ 

=  TV  ■ 


_N^r  1  Pa,rB o 

+  v'ff  RTr^r  ' 

UK,xJ  UKaJ  rr 


X(CpaWaN  a) 

Va 


-XlhaV-(WaNa)]-V-  q  racl  — V 

Va 


(64) 


where  Na-r  is  given  by  Faraday's  law.  Similarly,  in  Eq.  (21)  ®=0 
and  the  reaction  heat  is  introduced  as  boundary  heat  fluxes.  Thus,  Qrev.c  =  0 


'(  Q  revT  Q  irrev 1 

where 


(65) 


(66) 
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=  [AHr,(T  =  o/C)  +  nFVoc]  I  N  fuelra  ’  n  ra  1 

(67) 

c  =  ~  irciflactx  ^Iconx) 

(68) 

a  =  ~  ^rai^act^a^^lcon^a) 

(69) 

qrev  being  non-zero  only  at  the  fuel  electrode  reaction  surface,  it  is 
exothermic  if  SOFC  and  endothermic  if  SOEC;  and  q!rrev  being  non¬ 
zero  at  both  electrode  reaction  walls,  it  is  exothermic  at  both 
electrodes  and  for  both  operation  modes  (SORFC).  These  fluxes  are 
defined  normal  to  the  reaction  surfaces. 

Next,  the  channel  model  is  solved  for  both  the  fuel  and  air 
channels  using  the  SIMPLE  velocity-pressure-coupling  algorithm 
to  solve  the  mass-  and  momentum-conservation  equations.  The 
energy-conservation  equation  is,  at  this  point,  conveniently  solved 
for  the  entire  SORFC  domain  using  a  coupled  solver.  Information  is 
then  exchanged  at  the  channel-electrode  coupled  boundaries. 

Finally,  the  electrochemical  model  equations  are  solved  to  calcu¬ 
late  the  concentration  overpotentials  explicitly,  and  the  activation 
overpotentials  implicitly  using  a  Newton-Raphson  method.  Then,  the 
value  of  the  electric  potential  in  the  ion-conducting  phase  (electro¬ 
lyte)  at  the  reaction  boundaries  is  calculated  as  follows: 


RT 

SOFC:  <rc  =  2p  In 


Po2, 


Pref 


1/2' 


Vact,c 


(70) 


SOFC:  <ra=—  £»-_ln(^-  l+^  +  V 

H20,ra  ) 


(71) 


RT 

SOEC:  rp*rc=-—  In 
zt 


Pref 

Po2,rc 


1/2' 


+  Oact.c 


(72) 


SOEC:  cp*ra=  -E°  +  ^ln 

Zt 


PH20,r 


Ph2, 


tfact.a  ^ 


(73) 


These  values  are  set  as  at  the  corresponding  boundaries  of  the 
effective  electric-potential  field  in  the  electrolyte:  they  are  needed  for 
the  solution  of  the  charge-transfer  conservation  equation  in  the 
electrolyte  model,  which  is  solved  next,  thus  providing  the  current 
density  distribution  in  the  cell  at  the  preset  operating  voltage.  Mass 
and  heat  balances  are  finally  checked  to  ensure  conservation  at 
convergence. 


4.  Experimental  data 

The  applicability  and  validity  of  the  model  and  numerical  tool 
is  shown  in  the  following  sections  by  comparing  the  numerical 
results  with  the  data  measured  by  Laguna-Bercero  et  al.  [10]  for 
their  single,  self-supported,  hydrogen-fed,  microtubular  SORFC. 

The  experimental  test  rig  is  sketched  in  Fig.  3.  It  consists  of  a 
cylindrical  furnace,  at  the  centre  of  which  a  single  microtubular 
SOFC  is  placed.  The  furnace  has  a  metal  casing  and  a  ceramic 
cylindrical  wall  with  heating  resistances.  The  space  between  both 
surfaces  is  filled  with  a  thermally  insulating  material.  The  micro¬ 
tubular  SOFC  is  longer  than  the  furnace,  so  that  both  ends  are  kept 
away  from  the  heated  area  and  economical  and  simple  sealing 
mechanisms  can  be  used.  The  air  electrode,  which  is  shorter  than 
the  electrolyte/fuel-electrode  tube,  is  located  in  the  middle  of  the 
heated  perimeter  to  keep  the  temperature  of  the  cell  active  area 
close  as  possible  to  the  furnace  preset  one.  The  relevant  features  of 
the  experimental  test-rig  and  its  operating  conditions  are  sum¬ 
marised  in  Table  4. 

The  experimental  i-V  curves  reported  by  Laguna-Bercero  et  al. 
[10]  are  plotted  in  Fig.  4.  We  will  use  in  this  paper  their 


Metalic  casing 


Table  4 

Reference  microtubular  SORFC:  geometry  and  operating  conditions. 


Microtubular  geometry 

SORFC 

Fuel  channel  internal  radius  [10] 

1.2  mm 

Fuel  channel  length  [10] 

10  cm 

Anode  thickness  [10] 

400  pm 

Anode  length  [10] 

10  cm 

Anode  porosity  [54] 

£a 

41.4% 

Anode  tortuosity  factor  [64] 

3 

Anode  mean  pore  diameter  [54] 

dp, a 

1.2  pm 

Anode  permeability  [64] 

Bo, a 

7.65e-15  m2 

Electrolyte  thickness  [10] 

20  pm 

Electrolyte  length  [10] 

10  cm 

Cathode  thickness  [10] 

80  pm 

Cathode  length  [10] 

1  cm 

Cathode  porosity  [54] 

ec 

37.1% 

Cathode  tortuosity  factor  [64] 

*c 

5 

Cathode  mean  pore  diameter  [54] 

dp.c 

0.12  pm 

Cathode  permeability5 

Bo,c 

5.5e— 17  m2 

Air  channel  external  radiusa 

1  cm 

Air  channel  length5 

7.5  cm 

Solid  materials  properties 

SORFC 

Anode  thermal  conductivity  [65] 

x s.n 

3Wm4 K 1 

Electrolyte  thermal  conductivity  [65] 

2  W  m~'  K  1 

Electrolyte  emissivity3 

Sc 

0.2 

Ionic  conductivity  constant  [65] 

A, 

85  000  S  m_1 

Ionic  conductivity  constant  [65] 

Be 

11  000  K 

Cathode  thermal  conductivity  [65] 

2S,c 

4Wm-’  K-’ 

Cathode  emissivity  [66] 

Sc 

0.4 

Furnace  emissivityc 

Covert 

0.8 

Reference  operating  conditions 

SOFC 

SOEC 

Furnace  set  temperature  [10] 

Toven 

1168  K 

1168  K 

Operating  pressure  [10] 

Pout 

100  000  Pa 

100  000  Pa 

Fuel  temperature  at  inlet3 

T injiiel 

473  K 

473  K 

Fuel  composition,  H2  :  H20  :  N2[10]a 

A a,  in  file  1 

15:70:15 

15:70:15 

80:20:0 

16:20:64 

Fuel  flow-rate  (at  70%H20,  Troom)  [10] 

Qinjuel 

0.2  1  min-1 

0.2  1  min-1 

Fuel  flow-rate  (at  20%H20,Troom)a 

Qinfuel 

0.125  1  min-1 

0.125  1  min-1 

Air  composition,  02  :  N2 

Xa,in 

0.21:0.79 

0.21:0.79 

Air  flow- rate  (at  rroom)a 

Qin,air 

0.1 1  min~’ 

0.1 1  min-1 

a  Personal  communication. 

b  Estimated  with  the  Poiseuille-flow  equation  [67], 
c  Typical  value  assumed. 


experimental  results  for  the  SOEC  mode  up  to  1.2  V  since  the 
authors  subsequently  detected  electrolyte  degradation  for  higher 
voltages  [49],  Also,  the  experimental  open  circuit  voltage  for  the 
SOEC  at  20%  H20(  ~  0.87  V)  is  lower  that  what  it  is  expected 
according  to  the  Nernst  equation  ( ^  0.89  V).  Laguna-Bercero  et  al. 
confirm  that  this  mismatch  is  attributable  to  an  inaccurate 


1  Personal  communication. 
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SOFC-SOEC,  T  =  1 1 68  K 


Fig.  4.  Experimental  characteristic  curves  of  the  microtubular  SORFC  of  Laguna-Bercero  et  al.  [  10]. 


performance  of  the  humidifier.  The  actual  composition  of  the 
feeding  gases  has  been  instead  estimated  to  be  14%  H2,  30%  H20 
and  56%  N2.  This  is  therefore  the  composition  used  to  calculate  the 
numerical  curve  labelled  “SOECsim  :  pH20  =  0.2  atm". 


5.  Fitting  process 

In  this  section  we  describe  the  fitting  process  employed  to 
assign  values  to  the  electrochemical  parameters  of  a  particular  cell 
prior  to  using  the  model  for  its  simulation. 

The  need  of  this  fitting  process  stems  from  the  multiscale 
nature  of  the  several  phenomena  in  the  SORFC.  Thus  CFD  methods 
usually  resolve  the  macroscale  phenomena,  while  the  microscale 
processes  are  modelled  and  included  as  source  terms  in  the 
transport  equations.  This  is  the  case  of  the  electrochemical  reac¬ 
tion  rate,  given  by  the  kinetic  equations  of  Butler-Volmer, 
Eqs.  (50)  and  (51).  These  involve  the  following  electrochemical 
parameters:  ab, a.  «f,a.  “b.c.  «f,c.  io.a,  W,  the  values  of  which  are 
difficult  to  acquire  experimentally. 

In  the  CFD  modelling  literature,  the  charge  transfer  coefficients 
are  often  set  to  0.5  both  for  SOFC  [15,50,12]  and  SOEC  [17,19,51], 
This  corresponds  to  a  simple,  single-electron-exchange  reaction- 
mechanism  [46],  despite  the  fact  that  hydrogen-oxidation  reaction 
involves  2  electrons.  Nevertheless,  this  simplification  is  generally 
accepted  because  of  its  ease  of  use  and  apparently  good  results. 
The  present  authors  did  successfully  use  this  assumption  in  a 
previous  modelling  study  for  SOFCs  [16],  However,  we  find  that 
this  value  may  be  inappropriate  when  modelling  the  SOFC  and 
SOEC  modes  simultaneously,  as  shown  below. 

The  estimation  of  the  exchange  current  density  is  typically  the 
other  key  issue  in  the  parameter  fitting  process.  In  the  simplest 
approach  among  those  used  in  the  literature,  the  exchange  current 
density  is  assumed  constant  and  its  value  is  fitted,  as  in  [20],  In  more 
comprehensive  models,  experimental  correlations  are  used  to  account 
for  the  exchange-current-density  dependence  on  temperature  and 
species  concentrations.  Eqs.  (52)  and  (53)  represent  typical  correla¬ 
tions  for  the  exchange  current  densities  of  both  electrodes,  where  ya,  yc 


are  the  fitted  parameters  and  the  values  a,  b,  c  significantly  vary  from 
one  correlation  to  another  due  to  their  strong  dependence  on  the  cell 
materials  and  microstructure.  None  of  the  CFD  models  reviewed  by 
the  authors  using  one  of  these  empirical  correlations  from  the 
literature  justify  the  accuracy  of  the  chosen  correlation. 

In  the  following  subsections,  we  present  the  parameter  fitting 
process  performed  in  this  work  to  estimate  the  values  of  the 
charge-transfer  coefficients  and  the  exchange  current  densities. 
The  aim  is  to  provide  more  reliable  fitted  values  for  SOFC  and  SOEC 
that  can  be  used  in  a  quantitative  modelling  tool. 

5.3.  Estimation  of  the  charge  transfer  coefficients 

The  charge  transfer  coefficients  ( aba ,  af  a,  abc,  afx)  are  usually  set 
to  0.5  ( i.e .  single  step,  single  electron  transfer)  [15,50,12,17,19,51], 
for  the  reasons,  and  despite  the  shortcomings,  referred  to  above; 
therefore,  for  most  models  the  fitted  parameters  are  just  the 
exchange  current  densities. 

Because  of  our  previous  successful  experience  with  these 
values  for  SOFCs  [16],  we  have  in  the  first  instance  attempted  to 
use  them  for  the  present  SOFC-SOEC  model.  Fig.  5  shows  the 
experimental  data  of  Laguna-Bercero  et  al.  [10]  together  with  the 
numerical  results  obtained  when  taking  aba  =  aba  =  abx  =  afx  =  0.5 
(see  “case-0”  in  Table  5);  the  numerical  value  of  i  is  calculated  over 
the  air-electrode  active  area  as 

i  =  -T-[  ~t rc  ■  ~n  rc  dArc  (74) 

nrcjArc 

From  the  figure,  it  is  apparent  that  taking  all  the  charge  transfer 
coefficients  as  0.5  seems  to  work  for  the  SOFC  mode  (as  in  [16]) 
but  fails  to  model  the  SOEC  operation,  as  indicated  by  the  arrows 
in  magenta.  According  to  the  Butler-Volmer  equation  (Eq.  (50)), 
the  magnitude  of  the  current  density  depends  on  the  forward  and 
backward  charge  transfer  coefficients  (ab,a,af,a,«b.c,«f,c)  regardless 
of  whether  the  operating  mode  is  SOFC  or  SOEC  (forward  or 
backward).  This  is  not  satisfied  by  the  numerical  results  shown  in 
Fig.  5.  Therefore,  in  the  present  study,  we  reject  the  hypothesis 
ab  a  =  afa  =  abc  =  af  c  =  0.5  and  we  estimate  the  value  of  each 
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SOFC-SOEC,  T  =  1168K:aba  =  Ofa  =  abc  =  tXf  c  =  0.5 


SOECpxn:  70%H2O 
SOECoim  “  70%H2O 

■  70%H,O 


mm? 


exp- 


O 


SOFC, 

c 

SOEQ 


;e-0 
exp 
sim, case-0 


70%H2O 

20%H2O 

20%H2O 


SOFCp„n:  20%FI2O  O 
SOFCsimiCase.g:  20%H2O 


Fig.  5.  Characteristic  curves  of  a  microtubular  SORFC.  Experimental  data  (open  symbols)  for  the  cell  in  Table  4  [10]  and  simulation  results  (solid  symbols)  for  different  values 
of  ab,ii  =  «f,0  =  a-b,c  =  «f,c  (see  Table  5).  Arrows  highlight  discrepancies  between  measurements  and  modelling.  (For  interpretation  of  the  references  to  color  in  this  figure 
caption,  the  reader  is  referred  to  the  web  version  of  this  paper.) 


Table  5 

Electrochemical  parameters  for  both  operating  modes  and  for  the  different  cases  studied. 


Parameter 

SOFC 

Ref. 

SOEC 

Ref. 

SOFC 

Case-0 

SOEC 

Case-0 

SOFC 

Case-1 

SOEC 

Case-1 

SOFC 

Case-2 

SOEC 

Case-2 

SOFC 

Case-3 

SOEC 

Case-3 

SOFC 

Case-4 

SOEC 

Case-4 

<*b,a 

1 

0.1 

0.5 

0.5 

1 

0.1 

1 

0.1 

1 

0.1 

1 

0.1 

«f,a 

0.1 

1 

0.5 

0.5 

0.1 

1 

0.1 

1 

0.1 

1 

0.1 

1 

Eact.a 

120 

120 

120 

120 

120 

120 

120 

120 

120 

120 

120 

120 

Ya 

2.5e9 

2.5e9 

2.5e9 

2.5e9 

2.5e9 

2.5e9 

lelO 

lelO 

2.5e9 

2.5e9 

5e9 

5e9 

a 

0.11 

0.11 

0.11 

0.11 

1 

1 

1 

1 

1 

1 

1 

1 

b 

0.67 

0.67 

0.67 

0.67 

-0.5 

-0.5 

-0.5 

-0.5 

-0.5 

-0.5 

-0.5 

-0.5 

c 

0.25 

0.25 

0.25 

0.25 

0.25 

0.25 

0.25 

0.25 

0.25 

0.25 

0.25 

0.25 

ab,c 

0.4 

0.1 

0.5 

0.5 

0.4 

0.1 

0.4 

0.1 

0.4 

0.1 

0.4 

0.1 

«f,c 

0.1 

0.4 

0.5 

0.5 

0.1 

0.4 

0.1 

0.4 

0.1 

0.4 

0.1 

0.4 

Eact.c 

120 

120 

120 

120 

120 

120 

120 

120 

120 

120 

120 

120 

Yc 

2.6e9 

2.6e9 

2.6e9 

2.6e9 

2.6e9 

2.6e9 

2.6e9 

2.6e9 

lelO 

lelO 

1.8e9 

1.8e9 

charge  transfer  coefficient  that  best  fits  the  experimental  data.  The 
resulting  values  are  reported  in  Table  6;  for  SOFC  mode  these 
are:  ab-a  =  l,  «fa  =  0.1,  abc  =  0.4,  afc  =  0.1.  For  SOEC  mode  the 
forward  and  backward  directions  are  inverted  and  hence:  aba  = 
0.1,  at a  =  1,  abc  =  0.1,  «f  C  =  0.4. 

These  results  may  be  quantitatively  interpreted  according  to 
the  theory  of  the  Butler-Volmer  equation  for  multistep  reactions 
[46],  which  states  that 

«b+«t  =  -  (75) 

v 

with  n  being  the  number  of  electrons  exchanged  in  the  electro¬ 
chemical  reaction  and  v  the  number  of  occurrences  of  the  rate¬ 
determining  step  in  the  overall  reaction.  For  a  hydrogen  fed  SORFC, 
such  the  one  studied  in  this  work,  n=2.  Hence,  from  Eq.  (75)  vc  =  4 
and  va  =  1  -8,  meaning  that  the  limiting  step  in  the  air  electrode 
occurs  four  times  and  the  one  in  the  fuel  electrode  happens  almost 
twice  per  each  overall  reaction. 

These  values  agree  qualitatively  with  those  reported  by 
Zhu  and  Kee  [45]  (SOFC  mode:  aba  =  \.5,  afa  =  0.5,  ab.c  =  0.75, 
arc  =  0.25)  for  a  button  SOFC  with  similar  materials  as  those  used 


Table  6 

Reference  microtubular  SORFC:  electrochemical  parameters  for  both  operating 
modes. 


Electrochemical  parameters 

SOFC 

SOEC 

Fuel  transfer  coefficients3 

ab,a>  af,a 

1,  0.1 

0.1, 1 

Fuel  activation  energy  [48] 

Eact,a 

120  kj  mol-1 

120  kj  mol-1 

Fuel  pre-exponential  coefficient3 

Ya 

2.5e9Am-2 

2.5e9Am-2 

Exponent [47] 

a 

0.11 

0.11 

Exponent [47] 

b 

0.67 

0.67 

Exponent [48] 

c 

0.25 

0.25 

Air  transfer  coefficients3 

ab,c>  at jc 

0.4,  0.1 

0.1,  0.4 

Air  activation  energy  [48] 

Eact,c 

120  kj  mol-1 

120  kj  mol-1 

Air  pre-exponential  coefficient3 

Yc 

2.6e9  A  m-2 

2.6e9  A  m-2 

a  Fitted  parameters. 


by  Laguna-Bercero  et  al.  [10],  i.e.  Ni-YSZ,  YSZ,  LSM.  Our  results  also 
agree  with  the  experimental  data  of  Kawada  et  al.  [23]  for  a  Ni-YSZ 
SOFC  anode  (SOFC  mode:  aba  =  2,  af  a  =  1)  and  Kenney  and  Karan 
[24]  for  a  LSM/YSZ  SOFC  cathode  (ab,c  =  1.5,  afc  =  0.5).  The  agree¬ 
ment  is  particularly  revealing  in  two  respects:  (i)  the  forward 
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charge  transfer  coefficient  is  smaller  than  the  backward  one  for 
both  electrodes;  and  (ii)  the  charge  transfer  coefficients  are  larger 
in  the  fuel  electrode  than  in  the  air  one. 


A  corollary  of  this  research  is  that  the  often-employed  assump¬ 
tion  «b,a  =  «f,a  =  «b,c  =  ,c  =  0.5  fails  when  properly  validated  with 
a  complete  set  of  experimental  data  for  the  cell  operating  in  direct 
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Fig.  6.  Characteristic  curves  of  a  microtubular  SORFC.  Experimental  data  (open  symbols)  for  the  cell  in  Table  4  [10]  and  simulation  results  (solid  symbols,  see  Table  5): 
(a)  varying  i0,a;  (b)  varying  i0jC;  (c)  varying  simultaneously  i0>a  and  i0>c.  (For  interpretation  of  the  references  to  color  in  this  figure  caption,  the  reader  is  referred  to  the  web 
version  of  this  paper.) 
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and  reverse  mode  (SORFC).  If  the  model  is  only  validated  with 
SOFC  or  SOEC  data,  the  error  made  by  this  assumption  will  be 
compensated  by  the  data  fitting,  i.e.  the  exchange  current  densi¬ 
ties,  and  this  in  turn  can  over/under-estimate  the  activation 
overpotentials. 


5.2.  Estimation  of  the  exchange  current  density 

The  exchange  current  density  is  often  estimated  using  empirical 
correlations  such  Eqs.  (52)  and  (53),  where  a,  b,  c  are  set  by  the 
correlation  and  then  ya,  yc  are  the  only  fitted  parameters.  In  the  works 
reviewed,  c  is  generally  set  to  0.25  while  a  and  b  vary  significantly;  for 
instance  Yamamura's  model  specifies  a=l  and  b=  -  0.5  (used  in 
[48,28,52]),  Mogensen's  model  (used  in  [48,51,17])  chooses  a=b=l 
and  recent  experimental  works  report  lower  values  for  a  and  inter¬ 
mediate  values  for  b  (a=0.11,  b= 0.67)  [47,53], 

Yamamura's  values  were  used  in  our  previous  work  [16], 
Validation  was  restricted  to  a  single  set  of  SOFC  data,  and  an 
acceptable  agreement  was  found.  Thus  these  values  were  used  in  a 
first  attempt  in  the  present  work  to  simulate  a  similar  reversible 
cell  for  two  different  fuel  compositions.  Given  a= 1,  b=- 0.5 
(Yamamura)  and  c=0.25,  the  value  of  ya  and  yc  is  guessed  in  a 
trial  and  error  method  aiming  at  the  best  possible  agreement 
between  numerical  and  experimental  i-V  curves.  Representative 
results  from  this  fitting  process  are  plotted  in  Fig.  6. 

Fig.  6(a)  shows  the  experimental  data  [10]  (open  symbols)  and 
two  iterations  of  the  fitting  process,  labelled  “case-1”  and  “case-2” 
and  specified  in  Table  5.  As  indicated  in  the  figure  by  magenta 
arrows,  “case-1”  (solid  line  with  symbols)  underpredicts  the  mean 
current  density  of  the  SORFC  with  70%  H20  and  the  SOEC  with  20% 
H20;  on  the  contrary,  it  overpredicts  the  performance  when 
running  the  cell  as  an  SOFC  with  20%  H20.  In  order  to  get  a  better 
fitting,  another  value  of  ya  is  assumed  in  “case-2"  (dashed  line  with 
symbols).  The  results  of  “case-2”  improve  those  of  “case-1” 
reducing  the  performance  underestimation  of  the  latter  and 
providing  a  good  fit  for  the  whole  SORFC  operating  range  under 
conditions  of  70%  H20.  However,  as  it  is  pointed  out  with  dashed 
arrows  in  cyan,  “case-2”  overpredicts  the  current  density  for  both 
SOFC  and  SOEC  operation  with  20%  H20. 


Fig.  6(a)  shows  that  when  ya  is  modified  to  improve  the  fit  for 
curves  corresponding  to  high  steam  concentrations  in  the  fuel, 
then  the  fit  of  the  curves  at  low  steam  concentrations  worsens. 

Similar  to  Fig.  6(a),  Fig.  6(b)  shows  the  results  of  the  fitting  process 
as  the  value  of  yc  is  changed;  this  is  increased  from  “case-1"  (solid  line 
with  symbols)  to  “case-3”  (dashed  line  with  symbols).  It  is  found  that 
“case-3 "(see  cyan  dashed  arrows):  (i)  provides  a  good  fit  for  the  70%- 
H20  SOFC  and  the  20%-H20  SOEC  curves;  (ii)  slightly  improves  the 
agreement  between  the  70%-H2O  SOEC  curves;  and  (iii)  worsens  the 
agreement  of  the  20%-H20  SOFC  curves. 

Fig.  6(b)  shows  that  varying  the  value  of  yc  is  not  enough  to  get 
simultaneously  a  good  fit  between  numerical  and  experimental 
data  for  the  four  curves  considered. 

Finally,  Fig.  6(c)  shows  the  results  when  ya  and  yc  are  modified 
simultaneously.  Here  "case-4”  (dashed  line  with  symbols)  is  obtained 
with  a  lower  value  of  yc  and  a  higher  value  of  ya  than  those  used  in 
“case-1"  (solid  line  with  symbols).  The  20%-H2O  SORFC  performance  is 
well  predicted  in  “case-4"  but  it  fails  in  the  prediction  of  the  70  %-H20 
SORFC,  as  indicated  by  the  dashed  cyan  arrows. 

After  several  iterations  like  those  represented  by  “case-1”,  “case-2”, 
“case-3”  and  “case-4”,  we  conclude  that  there  is  no  pair  of  values  for  ya 
and  yc  that  makes  the  numerical  and  experimental  data  agree  simu¬ 
ltaneously  in  all  four  situations  studied.  However,  some  combinations 
provide  a  good  fit  of  the  SOFC  and  SOEC  curves  under  the  same 
operating  conditions,  e.g.  “case-2”  reproduces  the  performance  of  the 
70  %-H20  SORFC  and  “case-4”  that  of  the  20%-H2O  SORFC.  Therefore, 
we  conclude  that  the  model  can  describe  the  SORFC  performance  but 
not  its  dependence  on  the  fuel  concentration.  This  conclusion  points  to 
Yamamura's  values  as  the  cause  of  this  underperformance. 

The  work  of  Bieberle  et  al.  [53],  recently  supported  by  Utz  et  al. 
[47],  was  then  used  to  replace  Yamamura's  correlation.  A  fitting 
process  similar  to  that  presented  above  was  conducted  using 
Bieberle's  values  of  a  =  0.11  and  b= 0.67,  resulting  in  a  pair  of 
values  for  ya  and  yc  that  provide  good  fit  of  the  experimental  data 
both  for  low  and  high  steam  concentrations  and  in  both  SOFC  and 
SOEC  operations.  The  best  fit  is  that  shown  in  Fig.  7  for  the 
electrochemical  values  reported  in  Table  6. 

To  summarise  this  section,  we  have  shown  that  the  results  of  a 
macroscale  CFD  tool  strongly  depend  on  the  correlation  chosen  to 
model  the  dependence  of  the  exchange  current  density  with  the  gas 
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Fig.  7.  Characteristic  curves  of  a  microtubular  SORFC.  Experimental  data  (open  symbols)  [10]  and  simulation  results  (solid  symbols)  for  the  cell  in  Tables  4  and  6. 
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SOFC  @(H2:H20:N2  =  1 5:70:1 5  %,  T  =  1 1 68K) 
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Fig.  9.  Mean  exchange  current  density  over  the  electrodes  active  area  for  both  electrodes  (ia  a,  i0,c )  and  for  all  the  simulations  plotted  in  Fig.  7,  together  with  the  cell  average 
current  exchange  density  (i0,ceii  =  (io,a +ia,c)/2)  and  a  reference  experimental  value  [55], 


composition  at  the  microscale.  The  authors  recommend  the  use  of 
(at  least)  a  couple  of  experimental  data  sets  with  different  reactant 
compositions  to  check  whether  a  correlation  from  the  literature, 
expressing  this  dependence,  truly  applies  to  a  cell  being  modelled. 

6.  Validation 

The  validation  of  the  mathematical  model  and  numerical  tool 
presented  in  this  work  is  shown  in  Fig.  7,  where  the  experimental 


data  for  a  single,  self-supported,  hydrogen-fed,  microtubular  solid 
oxide  regenerative  fuel  cell  [10]  is  compared  with  the  numerical 
results  obtained  after  the  fitting  process  presented  in  the  former 
section.  The  values  of  the  fitted  parameters  are  reported  in  Table  6, 
namely:  (i)  the  pre-exponential  coefficients  (/ a ,  yc),  the  exponents 
(a,  fa)  in  Eq.  (52),  and  the  backward  and  forward  transfer  coeffi¬ 
cients  (ab  a,  ttf.o,  «b.c,  «f,c)- 

As  shown  in  Fig.  7,  numerical  and  experimental  curves  are  in 
good  agreement.  However,  the  applicability  of  such  numerical 
tool  must  be  further  verified  to  check  the  reliability  of  the  values 
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estimated  for  the  parameters.  In  the  literature,  however,  this 
additional  validation  is  rarely  performed,  probably  because  of 
the  lack  of  complementary  experimental  data  to  validate  against 
(e.g.  Electrochemical  Impedance  Spectroscopy,  EIS,  analysis). 

Fig.  8  shows  the  distribution  of  the  voltage  losses  obtained 
numerically  for  the  given  values  of  the  fitted  parameters.  It  is 
found  that  the  main  voltage  drop  is  due  to  the  activation  losses,  as 
it  may  be  presumed  from  experimental  data  [54,49].  Moreover, 
Fig.  9  shows  the  distribution  of  the  calculated  mean  exchange 
current  densities,  which  are  in  the  same  order  of  magnitude  as  the 
value  estimated  from  experimental  EIS  data  (ioceu  «  0.5  A  cm-2) 
for  a  similar  cell  running  as  SOFC  under  slightly  different  operating 
conditions  (V=0.95  V,  T=1168K,  H2:  H20=4.85:3)  [55],  It  is  thus 
shown  that  the  values  of  the  fitted  parameters  not  only  provide  a 
good  accord  between  the  numerical  and  experimental  character¬ 
istic  curves,  but  also  reliable  results  of  the  internal  physics  of 
the  cell. 


7.  Conclusions 

In  this  work,  we  have  presented  a  comprehensive  numerical 
tool  for  the  simulation  of  a  solid  oxide  regenerative  fuel  cell.  The 
tool  models  the  following  SORFC  phenomena:  (a)  mass  transport 
through  channels  and  porous  solids;  (b)  momentum  transport 
through  channels  and  porous  solids;  (c)  species  transfer  through 
channels  and  porous  media  (convection,  ordinary  diffusion,  Knud- 
sen  diffusion);  (d)  heat  transfer  in  gases  and  solids  by  conduction, 
convection  and  surface-to-surface  radiation;  (e)  charge  transport 
through  an  impervious  solid;  and  (f)  reversible  electrochemistry. 

The  validation  of  such  numerical  tool  involves  the  fitting  of 
the  unknown  parameters  in  the  electrochemical  model,  namely: 
(i)  the  pre-exponential  coefficients  (ya,  yc),  the  exponents  (a,  fa)  in 
Eq.  (52),  and  the  backward  and  forward  transfer  coefficients 

ac/a,  <*b,c>  otf,c )  • 

A  literature  review  regarding  this  fitting  process  is  also  pre¬ 
sented.  It  leads  to  the  conclusion  that  this  process  often  over¬ 
looked  in  the  literature,  where  it  is  (if  at  all  mentioned)  not  clearly 
explained  or  sufficiently  justified.  Thus,  the  fitting  is  here 
described  in  detail.  Two  corollaries  emerges  in  this  research. 

First,  the  often-employed  assumption  aba  =  «f =  abc  =  af  c  =  0.5 
fails  when  properly  validated  with  a  complete  set  of  experimental  data 
for  the  cell  operating  in  direct  and  reverse  mode  (SORFC).  Our  results 
suggest  that  the  forward  charge  transfer  coefficient  is  smaller  than  the 
backward  one  for  both  electrodes;  and  that  the  charge  transfer 
coefficients  are  larger  in  the  fuel  electrode  than  in  the  air  one. 

Second,  the  relevance  of  exchange-current-density  modelling  is 
also  highlighted  in  this  study;  we  have  shown  that  an  accurate 
modelling  of  the  dependence  of  the  exchange  current  density  on 
the  reactant  concentration  is  paramount  for  a  reliable  prediction  of 
the  overall  cell  performance. 

The  validity  of  the  numerical  results  provided  after  the  fitting 
process  for  some  of  the  uncertain  parameters  has  been  shown  by 
comparison  with  experimental  data  for  an  anode-supported 
microtubular  SORFC.  Good  agreement  is  found  not  only  between 
the  numerical  and  experimental  i-V  curves,  but  also  between 
numerical  and  experimental  electrochemical  parameters,  such  as 
activation  overpotentials  and  exchange  current  densities,  which 
supports  the  accuracy  of  the  fitted  values. 
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Appendix  A.  Thermophysical  properties 

A.l.  Viscosity,  p 


The  semi-empirical  formula  of  Wilke  is  used  to  estimate  the 
viscosity  of  a  multicomponent  fluid  [68]: 


P  = 


z 

Va 


Xqftq 

Sv/l(^ i) 


(A.l) 


where 

.  [i  +(/g/^)1/2(vvwu1/4]2 

al>  (8  +  8  WJWpj1'2 


(A.2) 


The  Sutherland  model  accounts  for  the  temperature  depen¬ 
dence  of  the  viscosity  pa  of  each  species  a : 


Pa  =  Poa 


f '  0Ba  + 

T  +  Sycr 


3/2 


(A.3) 


where  the  Sutherland-law  parameters  (>;0,T0,,,S^)  are  tabulated  for 
the  most  common  gases  in  [68], 


A.2.  Binary  diffusion  coefficient,  Da/1 


The  binary-diffusion  coefficients  are  modelled  according  to  the 
empirical  correlation  given  by  Wilke  [69]: 


Dap  =  2.628  x  10 


-3 


3  Wa  +  W„ 
2WaWp 


alpClDP 


(A.4) 


where  the  binary  diffusion  coefficient  Va/i  is  given  in  square 
centimetres  per  second;  the  temperature,  T,  in  Kelvin;  the  mole¬ 
cular  weight  of  species  a,  Wa,  in  kilogram  per  kilomole;  the 
pressure,  p,  in  atmospheres.  The  average  collision  diameter,  in 
angstroms,  is  calculated  as 


Oa/l  = 


a  + 

2 


(A.5) 


where  the  value  of  aa  (or  af)  is  tabulated  for  the  most  common 
species  in  [70],  £2D  is  the  dimensionless  collision  integral  in  the 
Lennard-Jones  12-6  potential  model: 

^  _  1.06036  0.19300  1.03587  1.76474 

^D—  -j-0  15610  +e(0.47635T)“*“g(1.52996T)“*“g(3.89411T)  (A.b) 


Here, 


(A.7) 


where  kB  is  the  Boltzmann  constant  and  e„  is  the  characteristic 
Lennard-Jones  energy  of  species  a.  For  the  most  common  gases  the 
value  of  ea  is  tabulated  in  [70]. 

To  convert  the  units  of  the  binaiy  diffusion  coefficient  and  the 
pressure  into  the  SI  system  (m2  s-1,  Pa),  Eq.  (A.4)  is  modified  as 
follows: 


Da/ 1  = 


2.628  x 

10.1325 - 


1f)-3  /T3W„  +  W( 
1U  y  1  2WaWp 

0%PdP 


(A.8) 


where  ua/J  is  the  only  parameter  that  is  not  expressed  in  SI  system 
units,  as  it  remains  in  angstroms. 
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The  diffusion  coefficient  of  a  species  a  in  a  multicomponent 
mixture  can  be  calculated  as 


The  reference  temperature  is  for  convenience  set  to  zero, 
Tref  =  0  K.  Hence,  from  Eqs.  (A.17)  and  (A.15): 


Tfrm  — 


1  — X 


a 


(A.9) 


h  n  T-l.  2“t2^  3“t3^ ,  U5»r5 

ha  =  —[a,aT+^T  +^-T  +^-T  +^T 


( A.  1 8 ) 


xa  being  the  molar  fraction  of  species  a  and  being  the  binary 

diffusion  coefficient  defined  above.  A  7.  Ionic  conductivity,  a 


A.3.  Knudsen  diffusion  coefficient,  DKa 


The  Kinetic  Theory  of  Gases  defines  the  Knudsen  diffusion 
coefficient  as  [71,72] 


(A.10) 


where  VKa  is  given  in  SI  units  if  the  temperature  is  expressed  in 
Kelvin,  the  molecular  weight  in  kilograms  per  kilomole,  the  mean 
pore  diameter  (dp)  in  meters  and  the  ideal  gas  law  constant  in 
Joule  per  kilomole  per  Kelvin. 


A.4.  Thermal  conductivity,  X 


Similar  to  viscosity  (Section  A.l),  the  thermal  conductivity  of 
the  fluid  is  given  by  Wilke's  formula  [68]: 


Va 


XaXa 


(All) 


where 

a  _[1  +(W^)1/2(W/»/WJ1/4]2 
0/J  (8  +  8  WJW/,?!2 


(A.12) 


The  thermal  conductivity  of  each  species  a,  Xa,  is  then  given  by  the 
Sutherland  model: 


where  the  Sutherland-law  parameters  (X0,T0i,S2)  are  tabulated  for 
the  most  common  gases  in  [68], 


A.5.  Specific  heat  at  constant  pressure,  Cp 

The  specific  heat  at  constant  pressure  of  a  multicomponent  gas 
is  calculated  as 

Cp  =  Sy^Cpo  (A. 14) 

Va 

where  Cpat  the  specific  heat  at  constant  pressure  of  species  a,  is 
provided  by  the  JANAF  thermochemical  tables: 

Cpa  =  ^-(a  i„  +  a2J + a3aT2  +  a4aT3  +  a5(J4)  (A.l  5) 

The  parameters  aJa,  a2a,  a3a,  a4a,  a5a  are  the  JANAF  constants  [73]. 


A.6.  Sensible  enthalpy,  h 

The  sensible  enthalpy  of  a  multicomponent  fluid  is  given  by 

h=2(vA)  (A.16) 

Va 

where  ha  is  the  sensible  enthalpy  of  each  component: 

K=  [T  Cpa  dT  (A.17) 

•JTref 


The  electrode  is  the  only  ion-conductor  SOFC-component. 
Theory  and  experiments  suggest  for  ion-conducting  electrolyte- 
materials  the  following  equation  for  the  ionic  conductivity: 


where  the  constants  Ae  and  Be  are  tabulated  for  most  common 
SOFC  electrolytes  in  [65]. 
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